#include "libmri.h"

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
using std::flush;

void eigsens3d (complex *smap, float *emap, int *imsize,
                complex *calib, int *csize, int *ksize,
                float tau, int qcols, int qriters)
{
	complex *V;
	cplxf *VVH;
	int nv;

	int im1 = imsize[0], im2 = imsize[1], im3 = imsize[2], nch = imsize[3];
	int ks1 = ksize[0], ks2 = ksize[1], ks3 = ksize[2];

	VVH = new cplxf [(2*ks1-1)*(2*ks2-1)*(2*ks3-1)*nch*nch];

	double ts, ts0 = timestamp();
	cout << "eigencalib3d_cheev() ... " << flush;
	calib3d_cheev (V, nv, calib, csize, ksize, tau);

	cout << ((ts = timestamp())-ts0) << endl;
	ts0 = ts;
	cout << "kreduce3d() ... " << flush;
	
	kreduce3d (VVH, (cplxf*) V, ksize, nch, nv);

	cout << ((ts = timestamp())-ts0) << endl;
	ts0 = ts;
	cout << "sens3d_orthiter() ... " << flush;

	calib3d_orthiter (smap, emap, (complex*) VVH, imsize, ksize, qcols, qriters);

	cout << ((ts = timestamp())-ts0) << endl;

	delete [] VVH;
	delete [] V;
}


#ifdef _MAIN_

int main (int argc, char **argv)
{
	shm_array_t *shm_calib = 0;
	int *xsize = 0, nxsize = 0; 
	int *ksize = 0, nksize = 0;
	int qriters = 4;
	int qcols = 2;
	char *smap_key = 0, *emap_key = 0;
	float tau = .1;

	struct arg* args [] = {
		new arg_t<int> ("qriters",  qriters),
		new arg_t<int> ("qcols",  qcols),
		new list_arg_t<int> ("imsize", xsize, nxsize, "x"),
		new list_arg_t<int> ("ksize", ksize, nksize, "x"),
		new arg_t<char*> ("smap", smap_key),
		new arg_t<char*> ("emap", emap_key),
		new arg_t<shm_array_t*> ("calib", shm_calib),
		new arg_t<float> ("tau", tau),
		0
	};

	get_args (argc, argv, args);

	if (xsize and nxsize != 4){
		cerr << "Must specify imsize=<int>x<int>x<int>x<int>" << endl;
		return -1;
	}

	if (shm_calib == 0){
		if (get_shm_array("calib", &shm_calib)){
			perror ("calib");
			return -1;
		}
	}
	if (shm_calib->type != shm_array_cplxf){
		cerr << "calib must be single-precision complex!" << endl;
		return -1;
	}
	if (shm_calib->ndims != 4){
		cerr << "calib must be <int>x<int>x<int>x<int>" << endl;
		return -1;
	}
	complex *calib = (complex*) shm_calib->data.c;

	int csize[4] = {shm_calib->dims[0], shm_calib->dims[1],
	                shm_calib->dims[2], shm_calib->dims[3]};

	if (!xsize){
		xsize = new int[nxsize = 4];
		xsize[0] = 64;
		xsize[1] = 64;
		xsize[2] = 64;
		xsize[3] = csize[3];
	}

	cout << "Image Size = "
		<< xsize[0] << " x " << xsize[1] << " x " << xsize[2] << " x "
		<< xsize[3] << endl;

	if (ksize and nksize != 3){
		cerr << "Must specify ksize=<int>x<int>x<int>" << endl;
		return -1;
	}

	if (!ksize){
		ksize = new int[3];
		ksize[0] = ksize[1] = ksize[2] = 7;
	}

	cout << "Kernel Size = "
		<< ksize[0] << " x " << ksize[1] << " x " << ksize[2] << endl;


	if (!smap_key)
		smap_key = strdup ("smap");

	shm_array_t sz_smap, *shm_smap;
	sz_smap.type = shm_array_cplxf;
	for (int i = 0; i < (sz_smap.ndims = 5)-1; ++i)
		sz_smap.dims[i] = xsize[i];
	sz_smap.dims[4] = qcols;

	if (make_shm_array (smap_key, &sz_smap, &shm_smap)){
		perror(smap_key);
		return -1;
	}
	complex *smap = (complex*) shm_smap->data.c;

	if (!emap_key)
		emap_key = strdup ("emap");

	shm_array_t sz_emap, *shm_emap;
	sz_emap.type = shm_array_float;
	for (int i = 0; i < (sz_emap.ndims = 4)-1; ++i)
		sz_emap.dims[i] = xsize[i];
	sz_emap.dims[3] = qcols;

	if (make_shm_array (emap_key, &sz_emap, &shm_emap)){
		perror(emap_key);
		return -1;
	}
	float *emap = (float *) shm_emap->data.f;

	if (not fftwf_init_threads()){
		cerr << "fftwf_init_threads() failure!" << endl;
		return -1;
	}

	eigsens3d (smap, emap, xsize, calib, csize, ksize, tau, qcols, qriters);

	return 0;
}

// ifdef _MAIN_
#endif
